home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / zunik.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  15.5 KB  |  325 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((zeror 0.0)
  12.       (zeroi 0.0)
  13.       (coner 1.0)
  14.       (conei 0.0)
  15.       (con (make-array 2 :element-type 'double-float))
  16.       (c (make-array 120 :element-type 'double-float)))
  17.   (declare (type (simple-array double-float (120)) c)
  18.            (type (simple-array double-float (2)) con)
  19.            (type double-float conei coner zeroi zeror))
  20.   (f2cl-lib:fset (f2cl-lib:fref con (1) ((1 2))) 0.3989422804014327)
  21.   (f2cl-lib:fset (f2cl-lib:fref con (2) ((1 2))) 1.2533141373155003)
  22.   (f2cl-lib:fset (f2cl-lib:fref c (1) ((1 120))) 1.0)
  23.   (f2cl-lib:fset (f2cl-lib:fref c (2) ((1 120))) -0.20833333333333337)
  24.   (f2cl-lib:fset (f2cl-lib:fref c (3) ((1 120))) 0.125)
  25.   (f2cl-lib:fset (f2cl-lib:fref c (4) ((1 120))) 0.3342013888888889)
  26.   (f2cl-lib:fset (f2cl-lib:fref c (5) ((1 120))) -0.40104166666666674)
  27.   (f2cl-lib:fset (f2cl-lib:fref c (6) ((1 120))) 0.0703125)
  28.   (f2cl-lib:fset (f2cl-lib:fref c (7) ((1 120))) -1.0258125964506173)
  29.   (f2cl-lib:fset (f2cl-lib:fref c (8) ((1 120))) 1.8464626736111112)
  30.   (f2cl-lib:fset (f2cl-lib:fref c (9) ((1 120))) -0.8912109375)
  31.   (f2cl-lib:fset (f2cl-lib:fref c (10) ((1 120))) 0.0732421875)
  32.   (f2cl-lib:fset (f2cl-lib:fref c (11) ((1 120))) 4.669584423426247)
  33.   (f2cl-lib:fset (f2cl-lib:fref c (12) ((1 120))) -11.207002616222994)
  34.   (f2cl-lib:fset (f2cl-lib:fref c (13) ((1 120))) 8.78912353515625)
  35.   (f2cl-lib:fset (f2cl-lib:fref c (14) ((1 120))) -2.3640869140625)
  36.   (f2cl-lib:fset (f2cl-lib:fref c (15) ((1 120))) 0.112152099609375)
  37.   (f2cl-lib:fset (f2cl-lib:fref c (16) ((1 120))) -28.212072558200244)
  38.   (f2cl-lib:fset (f2cl-lib:fref c (17) ((1 120))) 84.63621767460074)
  39.   (f2cl-lib:fset (f2cl-lib:fref c (18) ((1 120))) -91.81824154324002)
  40.   (f2cl-lib:fset (f2cl-lib:fref c (19) ((1 120))) 42.53499874538846)
  41.   (f2cl-lib:fset (f2cl-lib:fref c (20) ((1 120))) -7.368794359479632)
  42.   (f2cl-lib:fset (f2cl-lib:fref c (21) ((1 120))) 0.22710800170898437)
  43.   (f2cl-lib:fset (f2cl-lib:fref c (22) ((1 120))) 212.57013003921713)
  44.   (f2cl-lib:fset (f2cl-lib:fref c (23) ((1 120))) -765.2524681411817)
  45.   (f2cl-lib:fset (f2cl-lib:fref c (24) ((1 120))) 1059.9904525279999)
  46.   (f2cl-lib:fset (f2cl-lib:fref c (25) ((1 120))) -699.5796273761325)
  47.   (f2cl-lib:fset (f2cl-lib:fref c (26) ((1 120))) 218.1905117442116)
  48.   (f2cl-lib:fset (f2cl-lib:fref c (27) ((1 120))) -26.491430486951554)
  49.   (f2cl-lib:fset (f2cl-lib:fref c (28) ((1 120))) 0.5725014209747314)
  50.   (f2cl-lib:fset (f2cl-lib:fref c (29) ((1 120))) -1919.457662318407)
  51.   (f2cl-lib:fset (f2cl-lib:fref c (30) ((1 120))) 8061.722181737309)
  52.   (f2cl-lib:fset (f2cl-lib:fref c (31) ((1 120))) -13586.550006434138)
  53.   (f2cl-lib:fset (f2cl-lib:fref c (32) ((1 120))) 11655.393336864532)
  54.   (f2cl-lib:fset (f2cl-lib:fref c (33) ((1 120))) -5305.646978613403)
  55.   (f2cl-lib:fset (f2cl-lib:fref c (34) ((1 120))) 1200.9029132163525)
  56.   (f2cl-lib:fset (f2cl-lib:fref c (35) ((1 120))) -108.09091978839464)
  57.   (f2cl-lib:fset (f2cl-lib:fref c (36) ((1 120))) 1.7277275025844574)
  58.   (f2cl-lib:fset (f2cl-lib:fref c (37) ((1 120))) 20204.291330966145)
  59.   (f2cl-lib:fset (f2cl-lib:fref c (38) ((1 120))) -96980.5983886375)
  60.   (f2cl-lib:fset (f2cl-lib:fref c (39) ((1 120))) 192547.00123253153)
  61.   (f2cl-lib:fset (f2cl-lib:fref c (40) ((1 120))) -203400.17728041552)
  62.   (f2cl-lib:fset (f2cl-lib:fref c (41) ((1 120))) 122200.46498301746)
  63.   (f2cl-lib:fset (f2cl-lib:fref c (42) ((1 120))) -41192.65496889756)
  64.   (f2cl-lib:fset (f2cl-lib:fref c (43) ((1 120))) 7109.514302489364)
  65.   (f2cl-lib:fset (f2cl-lib:fref c (44) ((1 120))) -493.915304773088)
  66.   (f2cl-lib:fset (f2cl-lib:fref c (45) ((1 120))) 6.074042001273483)
  67.   (f2cl-lib:fset (f2cl-lib:fref c (46) ((1 120))) -242919.1879005513)
  68.   (f2cl-lib:fset (f2cl-lib:fref c (47) ((1 120))) 1311763.6146629772)
  69.   (f2cl-lib:fset (f2cl-lib:fref c (48) ((1 120))) -2998015.9185381066)
  70.   (f2cl-lib:fset (f2cl-lib:fref c (49) ((1 120))) 3763271.297656404)
  71.   (f2cl-lib:fset (f2cl-lib:fref c (50) ((1 120))) -2813563.226586534)
  72.   (f2cl-lib:fset (f2cl-lib:fref c (51) ((1 120))) 1268365.2733216248)
  73.   (f2cl-lib:fset (f2cl-lib:fref c (52) ((1 120))) -331645.1724845636)
  74.   (f2cl-lib:fset (f2cl-lib:fref c (53) ((1 120))) 45218.76898136273)
  75.   (f2cl-lib:fset (f2cl-lib:fref c (54) ((1 120))) -2499.8304818112097)
  76.   (f2cl-lib:fset (f2cl-lib:fref c (55) ((1 120))) 24.380529699556064)
  77.   (f2cl-lib:fset (f2cl-lib:fref c (56) ((1 120))) 3284469.853072038)
  78.   (f2cl-lib:fset (f2cl-lib:fref c (57) ((1 120))) -1.9706819118432228e+7)
  79.   (f2cl-lib:fset (f2cl-lib:fref c (58) ((1 120))) 5.095260249266464e+7)
  80.   (f2cl-lib:fset (f2cl-lib:fref c (59) ((1 120))) -7.410514821153266e+7)
  81.   (f2cl-lib:fset (f2cl-lib:fref c (60) ((1 120))) 6.634451227472903e+7)
  82.   (f2cl-lib:fset (f2cl-lib:fref c (61) ((1 120))) -3.756717666076335e+7)
  83.   (f2cl-lib:fset (f2cl-lib:fref c (62) ((1 120))) 1.3288767166421817e+7)
  84.   (f2cl-lib:fset (f2cl-lib:fref c (63) ((1 120))) -2785618.1280864547)
  85.   (f2cl-lib:fset (f2cl-lib:fref c (64) ((1 120))) 308186.4046126624)
  86.   (f2cl-lib:fset (f2cl-lib:fref c (65) ((1 120))) -13886.089753717042)
  87.   (f2cl-lib:fset (f2cl-lib:fref c (66) ((1 120))) 110.01714026924674)
  88.   (f2cl-lib:fset (f2cl-lib:fref c (67) ((1 120))) -4.932925366450996e+7)
  89.   (f2cl-lib:fset (f2cl-lib:fref c (68) ((1 120))) 3.2557307418576575e+8)
  90.   (f2cl-lib:fset (f2cl-lib:fref c (69) ((1 120))) -9.394623596815784e+8)
  91.   (f2cl-lib:fset (f2cl-lib:fref c (70) ((1 120))) 1.55359689957058e+9)
  92.   (f2cl-lib:fset (f2cl-lib:fref c (71) ((1 120))) -1.6210805521083374e+9)
  93.   (f2cl-lib:fset (f2cl-lib:fref c (72) ((1 120))) 1.1068428168230145e+9)
  94.   (f2cl-lib:fset (f2cl-lib:fref c (73) ((1 120))) -4.958897842750303e+8)
  95.   (f2cl-lib:fset (f2cl-lib:fref c (74) ((1 120))) 1.4206290779753308e+8)
  96.   (f2cl-lib:fset (f2cl-lib:fref c (75) ((1 120))) -2.4474062725738724e+7)
  97.   (f2cl-lib:fset (f2cl-lib:fref c (76) ((1 120))) 2243768.1779224495)
  98.   (f2cl-lib:fset (f2cl-lib:fref c (77) ((1 120))) -84005.4336030241)
  99.   (f2cl-lib:fset (f2cl-lib:fref c (78) ((1 120))) 551.3358961220206)
  100.   (f2cl-lib:fset (f2cl-lib:fref c (79) ((1 120))) 8.147890961183122e+8)
  101.   (f2cl-lib:fset (f2cl-lib:fref c (80) ((1 120))) -5.866481492051847e+9)
  102.   (f2cl-lib:fset (f2cl-lib:fref c (81) ((1 120))) 1.8688207509295826e+10)
  103.   (f2cl-lib:fset (f2cl-lib:fref c (82) ((1 120))) -3.4632043388158773e+10)
  104.   (f2cl-lib:fset (f2cl-lib:fref c (83) ((1 120))) 4.1280185579753975e+10)
  105.   (f2cl-lib:fset (f2cl-lib:fref c (84) ((1 120))) -3.302659974980072e+10)
  106.   (f2cl-lib:fset (f2cl-lib:fref c (85) ((1 120))) 1.79542137311556e+10)
  107.   (f2cl-lib:fset (f2cl-lib:fref c (86) ((1 120))) -6.563293792619285e+9)
  108.   (f2cl-lib:fset (f2cl-lib:fref c (87) ((1 120))) 1.5592798648792575e+9)
  109.   (f2cl-lib:fset (f2cl-lib:fref c (88) ((1 120))) -2.251056618894153e+8)
  110.   (f2cl-lib:fset (f2cl-lib:fref c (89) ((1 120))) 1.7395107553978165e+7)
  111.   (f2cl-lib:fset (f2cl-lib:fref c (90) ((1 120))) -549842.3275722887)
  112.   (f2cl-lib:fset (f2cl-lib:fref c (91) ((1 120))) 3038.0905109223845)
  113.   (f2cl-lib:fset (f2cl-lib:fref c (92) ((1 120))) -1.4679261247695616e+10)
  114.   (f2cl-lib:fset (f2cl-lib:fref c (93) ((1 120))) 1.1449823773202579e+11)
  115.   (f2cl-lib:fset (f2cl-lib:fref c (94) ((1 120))) -3.990961752244665e+11)
  116.   (f2cl-lib:fset (f2cl-lib:fref c (95) ((1 120))) 8.192186695485775e+11)
  117.   (f2cl-lib:fset (f2cl-lib:fref c (96) ((1 120))) -1.0983751560812234e+12)
  118.   (f2cl-lib:fset (f2cl-lib:fref c (97) ((1 120))) 1.008158106865382e+12)
  119.   (f2cl-lib:fset (f2cl-lib:fref c (98) ((1 120))) -6.453648692453764e+11)
  120.   (f2cl-lib:fset (f2cl-lib:fref c (99) ((1 120))) 2.879006499061506e+11)
  121.   (f2cl-lib:fset (f2cl-lib:fref c (100) ((1 120))) -8.786707217802326e+10)
  122.   (f2cl-lib:fset (f2cl-lib:fref c (101) ((1 120))) 1.763473060683497e+10)
  123.   (f2cl-lib:fset (f2cl-lib:fref c (102) ((1 120))) -2.167164983223795e+9)
  124.   (f2cl-lib:fset (f2cl-lib:fref c (103) ((1 120))) 1.4315787671888897e+8)
  125.   (f2cl-lib:fset (f2cl-lib:fref c (104) ((1 120))) -3871833.4425726123)
  126.   (f2cl-lib:fset (f2cl-lib:fref c (105) ((1 120))) 18257.755474293175)
  127.   (f2cl-lib:fset (f2cl-lib:fref c (106) ((1 120))) 2.8646403571767903e+11)
  128.   (f2cl-lib:fset (f2cl-lib:fref c (107) ((1 120))) -2.406297900028504e+12)
  129.   (f2cl-lib:fset (f2cl-lib:fref c (108) ((1 120))) 9.109341185239899e+12)
  130.   (f2cl-lib:fset (f2cl-lib:fref c (109) ((1 120))) -2.051689941093444e+13)
  131.   (f2cl-lib:fset (f2cl-lib:fref c (110) ((1 120))) 3.056512551993532e+13)
  132.   (f2cl-lib:fset (f2cl-lib:fref c (111) ((1 120))) -3.166708858478516e+13)
  133.   (f2cl-lib:fset (f2cl-lib:fref c (112) ((1 120))) 2.334836404458184e+13)
  134.   (f2cl-lib:fset (f2cl-lib:fref c (113) ((1 120))) -1.2320491305598287e+13)
  135.   (f2cl-lib:fset (f2cl-lib:fref c (114) ((1 120))) 4.6127257808491323e+12)
  136.   (f2cl-lib:fset (f2cl-lib:fref c (115) ((1 120))) -1.1965528801961815e+12)
  137.   (f2cl-lib:fset (f2cl-lib:fref c (116) ((1 120))) 2.0591450323241003e+11)
  138.   (f2cl-lib:fset (f2cl-lib:fref c (117) ((1 120))) -2.182292775752922e+10)
  139.   (f2cl-lib:fset (f2cl-lib:fref c (118) ((1 120))) 1.2470092935127104e+9)
  140.   (f2cl-lib:fset (f2cl-lib:fref c (119) ((1 120))) -2.918838812222081e+7)
  141.   (f2cl-lib:fset (f2cl-lib:fref c (120) ((1 120))) 118838.42625678325)
  142.   (defun zunik
  143.          (zrr zri fnu ikflg ipmtr tol init phir phii zeta1r zeta1i zeta2r
  144.           zeta2i sumr sumi cwrkr cwrki)
  145.     (declare (type f2cl-lib:integer4 ikflg ipmtr init)
  146.              (type double-float zrr zri fnu tol phir phii zeta1r zeta1i zeta2r
  147.               zeta2i sumr sumi)
  148.              (type (simple-array double-float (*)) cwrkr cwrki))
  149.     (prog ((i 0) (idum 0) (j 0) (k 0) (l 0) (ac 0.0) (crfni 0.0) (crfnr 0.0)
  150.            (rfn 0.0) (si 0.0) (sr 0.0) (sri 0.0) (srr 0.0) (sti 0.0) (str 0.0)
  151.            (test 0.0) (ti 0.0) (tr 0.0) (t2i 0.0) (t2r 0.0) (zni 0.0)
  152.            (znr 0.0))
  153.       (declare
  154.        (type double-float znr zni t2r t2i tr ti test str sti srr sri sr si rfn
  155.         crfnr crfni ac)
  156.        (type f2cl-lib:integer4 l k j idum i))
  157.       (if (/= init 0) (go label40))
  158.       (setf rfn (/ 1.0 fnu))
  159.       (setf test (* (f2cl-lib:d1mach 1) 1000.0))
  160.       (setf ac (* fnu test))
  161.       (if (or (> (abs zrr) ac) (> (abs zri) ac)) (go label15))
  162.       (setf zeta1r (+ (* 2.0 (abs (f2cl-lib:flog test))) fnu))
  163.       (setf zeta1i 0.0)
  164.       (setf zeta2r fnu)
  165.       (setf zeta2i 0.0)
  166.       (setf phir 1.0)
  167.       (setf phii 0.0)
  168.       (go end_label)
  169.      label15
  170.       (setf tr (* zrr rfn))
  171.       (setf ti (* zri rfn))
  172.       (setf sr (+ coner (- (* tr tr) (* ti ti))))
  173.       (setf si (+ conei (+ (* tr ti) (* ti tr))))
  174.       (multiple-value-bind
  175.           (var-0 var-1 var-2 var-3)
  176.           (zsqrt sr si srr sri)
  177.         (declare (ignore var-0 var-1))
  178.         (setf srr var-2)
  179.         (setf sri var-3))
  180.       (setf str (+ coner srr))
  181.       (setf sti (+ conei sri))
  182.       (multiple-value-bind
  183.           (var-0 var-1 var-2 var-3 var-4 var-5)
  184.           (zdiv str sti tr ti znr zni)
  185.         (declare (ignore var-0 var-1 var-2 var-3))
  186.         (setf znr var-4)
  187.         (setf zni var-5))
  188.       (multiple-value-bind
  189.           (var-0 var-1 var-2 var-3 var-4)
  190.           (zlog znr zni str sti idum)
  191.         (declare (ignore var-0 var-1))
  192.         (setf str var-2)
  193.         (setf sti var-3)
  194.         (setf idum var-4))
  195.       (setf zeta1r (* fnu str))
  196.       (setf zeta1i (* fnu sti))
  197.       (setf zeta2r (* fnu srr))
  198.       (setf zeta2i (* fnu sri))
  199.       (multiple-value-bind
  200.           (var-0 var-1 var-2 var-3 var-4 var-5)
  201.           (zdiv coner conei srr sri tr ti)
  202.         (declare (ignore var-0 var-1 var-2 var-3))
  203.         (setf tr var-4)
  204.         (setf ti var-5))
  205.       (setf srr (* tr rfn))
  206.       (setf sri (* ti rfn))
  207.       (multiple-value-bind
  208.           (var-0 var-1 var-2 var-3)
  209.           (zsqrt srr sri (f2cl-lib:fref cwrkr (16) ((1 16)))
  210.            (f2cl-lib:fref cwrki (16) ((1 16))))
  211.         (declare (ignore var-0 var-1))
  212.         (f2cl-lib:fset (f2cl-lib:fref cwrkr (16) ((1 16))) var-2)
  213.         (f2cl-lib:fset (f2cl-lib:fref cwrki (16) ((1 16))) var-3))
  214.       (setf phir
  215.               (* (f2cl-lib:fref cwrkr (16) ((1 16)))
  216.                  (f2cl-lib:fref con (ikflg) ((1 2)))))
  217.       (setf phii
  218.               (* (f2cl-lib:fref cwrki (16) ((1 16)))
  219.                  (f2cl-lib:fref con (ikflg) ((1 2)))))
  220.       (if (/= ipmtr 0) (go end_label))
  221.       (multiple-value-bind
  222.           (var-0 var-1 var-2 var-3 var-4 var-5)
  223.           (zdiv coner conei sr si t2r t2i)
  224.         (declare (ignore var-0 var-1 var-2 var-3))
  225.         (setf t2r var-4)
  226.         (setf t2i var-5))
  227.       (f2cl-lib:fset (f2cl-lib:fref cwrkr (1) ((1 16))) coner)
  228.       (f2cl-lib:fset (f2cl-lib:fref cwrki (1) ((1 16))) conei)
  229.       (setf crfnr coner)
  230.       (setf crfni conei)
  231.       (setf ac 1.0)
  232.       (setf l 1)
  233.       (f2cl-lib:fdo (k 2 (f2cl-lib:int-add k 1))
  234.                     ((> k 15) nil)
  235.         (tagbody
  236.           (setf sr zeror)
  237.           (setf si zeroi)
  238.           (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
  239.                         ((> j k) nil)
  240.             (tagbody
  241.               (setf l (f2cl-lib:int-add l 1))
  242.               (setf str
  243.                       (+ (- (* sr t2r) (* si t2i))
  244.                          (f2cl-lib:fref c (l) ((1 120)))))
  245.               (setf si (+ (* sr t2i) (* si t2r)))
  246.               (setf sr str)
  247.              label10))
  248.           (setf str (- (* crfnr srr) (* crfni sri)))
  249.           (setf crfni (+ (* crfnr sri) (* crfni srr)))
  250.           (setf crfnr str)
  251.           (f2cl-lib:fset (f2cl-lib:fref cwrkr (k) ((1 16)))
  252.                          (- (* crfnr sr) (* crfni si)))
  253.           (f2cl-lib:fset (f2cl-lib:fref cwrki (k) ((1 16)))
  254.                          (+ (* crfnr si) (* crfni sr)))
  255.           (setf ac (* ac rfn))
  256.           (setf test
  257.                   (coerce
  258.                    (+ (abs (f2cl-lib:fref cwrkr (k) ((1 16))))
  259.                       (abs (f2cl-lib:fref cwrki (k) ((1 16)))))
  260.                    'double-float))
  261.           (if (and (< ac tol) (< test tol)) (go label30))
  262.          label20))
  263.       (setf k 15)
  264.      label30
  265.       (setf init k)
  266.      label40
  267.       (if (= ikflg 2) (go label60))
  268.       (setf sr zeror)
  269.       (setf si zeroi)
  270.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  271.                     ((> i init) nil)
  272.         (tagbody
  273.           (setf sr (+ sr (f2cl-lib:fref cwrkr (i) ((1 16)))))
  274.           (setf si (+ si (f2cl-lib:fref cwrki (i) ((1 16)))))
  275.          label50))
  276.       (setf sumr sr)
  277.       (setf sumi si)
  278.       (setf phir
  279.               (* (f2cl-lib:fref cwrkr (16) ((1 16)))
  280.                  (f2cl-lib:fref con (1) ((1 2)))))
  281.       (setf phii
  282.               (* (f2cl-lib:fref cwrki (16) ((1 16)))
  283.                  (f2cl-lib:fref con (1) ((1 2)))))
  284.       (go end_label)
  285.      label60
  286.       (setf sr zeror)
  287.       (setf si zeroi)
  288.       (setf tr coner)
  289.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  290.                     ((> i init) nil)
  291.         (tagbody
  292.           (setf sr (+ sr (* tr (f2cl-lib:fref cwrkr (i) ((1 16))))))
  293.           (setf si (+ si (* tr (f2cl-lib:fref cwrki (i) ((1 16))))))
  294.           (setf tr (- tr))
  295.          label70))
  296.       (setf sumr sr)
  297.       (setf sumi si)
  298.       (setf phir
  299.               (* (f2cl-lib:fref cwrkr (16) ((1 16)))
  300.                  (f2cl-lib:fref con (2) ((1 2)))))
  301.       (setf phii
  302.               (* (f2cl-lib:fref cwrki (16) ((1 16)))
  303.                  (f2cl-lib:fref con (2) ((1 2)))))
  304.       (go end_label)
  305.      end_label
  306.       (return
  307.        (values nil
  308.                nil
  309.                nil
  310.                nil
  311.                nil
  312.                nil
  313.                init
  314.                phir
  315.                phii
  316.                zeta1r
  317.                zeta1i
  318.                zeta2r
  319.                zeta2i
  320.                sumr
  321.                sumi
  322.                nil
  323.                nil)))))
  324.  
  325.